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ABSTRACT 

A viscous instability in shearing laminar axisymmetric hydrodynamic flows around a gravi¬ 
tating center is described. In the linearized hydrodynamic equations written in the Boussinesq 
approximation with microscopic molecular transport coefficients, the instability arises when 
the viscous dissipation is taken into account in the energy equation. Using the local WKB 
approximation, we derive a third-order algebraic dispersion equation with two modes rep¬ 
resenting the modified Rayleigh modes R-h and R-, and the third X-mode. We show that in 
thin accretion flows the viscosity destabilizes one of the Rayleigh modes in a wide range of 
wavenumbers, while the X-mode always remains stable. In Keplerian flows, the instability 
increment is found to be a few Keplerian rotational periods at wavelengths with kr ~ 10 - 50. 
This instability may cause turbulence in astrophysical accretion discs even in the absence of 
magnetic field. 
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1 INTRODUCTION 

The origin of turbulence in accretion discs is an outstanding prob¬ 
lem in astrophysics . The dimensionless phenom enological parame¬ 
ter a introduced bv ishakura & SunvaevI l ll973t) for assumed turbu¬ 
lent eddy viscosity and chaotic magnetic fields turned out to be very 
useful in describing physical properties of accretion discs. Analy¬ 
sis of different observations (e. g., the behaviour of non-stationary 
accre tion discs in X-ray novae jSuleimano v. Lipunova & Shakural 
l2008h and dwarf-nova and AM CVn stars jKotko & Laso32012l) ) 
suggest a rather large values a ~ 0.3, indicating the presence of 
well-developed turbulence in the disc. In Keplerian accretion discs, 
the angular momentum increases with radius, making the flow sta¬ 
ble against small hydrodynamic perturbations according to the clas¬ 
sical Rayleigh criterion. When the small magnetic field is present 
in fully ionized gas, a popular mechanism quenching the insta¬ 
bility is t he Velikhov-Chandrasekhar magneto-rotational instabil¬ 
ity (M RI) I Wli^o^[l95^ Chandrasekhaall9^ : lBalbus & HawlevI 
I l99lh (see lBalbus & Hawi^ h998l) for a detailed review). In spite 
of being a powerful instab ility, MRI has its own limitations (see 
e.g. lGoodman_&_Xjj[l99^, for di scussion of parasiting instabili¬ 
ties and Shakura & Postno^ ^20141) . for discussion of applications 
to thin Keplerian accretion discs). 

As for the purely hydrodynamic case, so far there has been no 
clear criterion of the hydrodynamic turbulence. In a Keplerian flow, 
there are different mechanisms for small perturbations growth, such 
as linear growth of transient pe rturbations (see the recent study 
Izhuravlev & RazdoburdinI ( |2014|) and references therein), but the 
transition of these perturbations to turbulence, which is a strongly 
non-linear process, remains unclear. 
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In this paper we perform the linear stability analysis of shear¬ 
ing laminar hydrodynamic viscous flows with arbitrary rotation 
laws in the form oc taking into account the viscous heating 
and thermal conductivity in the energy equation. Unlike many pre¬ 
vious works, we use microscopic molecular transport coetficients to 
describe the viscosity and heat conductivity. These terms in the en¬ 
ergy equation make one of the Rayleigh modes unstable (i.e. make 
their amplitude exponentially growing) in a wide range of wave 
numbers for perturbations normal to the direction of the wave vec¬ 
tor. Physically, the instability may be due to the viscously heated 
gas being unstable to convection in the gravity field of the central 
object in the absence of the background entropy gradients. 

The instability increment decreases (but does not vanish) with 
increasing thermal conductivity and is maximum in cold neutral 
flows with largest Prandtl numbers. We discuss the relevance of the 
found viscous instability to the generation of turbulence in laminar 
thin accretion flows. 

The pulsati onal in s tabilit y of viscous accretion discs was 
first studied by iKatd h978l) . The viscous instability of the 
standard turbulent Shakura-S unyaev a-discs was inv e stigate d 
in many papers (see , e.g., |Blumenthalj^in&Jft|nd ( Il984h : 
iKlev. Papaloizou & LinI h993h : lLatter & Qgilviel 1 20od) , among 


others). (Note that in the latter papers the viscous instability is re¬ 
ferred to as ’viscous overstability’, i.e. when arising as an exponen¬ 
tially gro wing oscillations, in analogy with stellar pulsations dis¬ 
cussed by lEddingtonl l ll926h . and the term ’instability’ is reserved 
for purely imaginary negative modes.) 

In Section[2]we derive the basic dispersion equation. To make 
the physical case as simple as possible, we work in the Boussinesq 
approximation (i.e. consider the fluid incompressible V« = 0, keep 
the Eulerian pressure variations non-zero only in the equations of 
motion and put them zero in the energy equation) and take small 
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perturbations in the form of plane waves in the direction transver¬ 
sal to the wave propagation. For such perturbations we derive a 
third-order algebaraic dispersion equation with imaginary coeffi¬ 
cients, which has three solutions: one Rayleigh mode with positive 
real part (R+), the Rayleigh mode with negative real part (R-) and a 
new mode which has zero real and imaginary part at ^ 0 (the X- 
mode). We find that one of the Rayleigh modes becomes unstable 
(exponentially growing) at large wavelengths, while the X-mode 
remains stable at all wavelengths. In Section[3]we analyze the ob¬ 
tained the dispersion equation. First we rewrite it in the dimension¬ 
less form, and then perform its numerical analysis for several im¬ 
portant cases of thermal conductivity (for purely electron conduc¬ 
tivity in the fully ionized plasma, for the case where the radiation 
conductivity is important, and for the case of neutral monoatomic 
hydrogen gas where the heat conductivity and viscosity are caused 
by the same particles). In Section[5]we discuss the applicability of 
the approximation of incompressibility we use and damping of the 
instability due to entropy gradients in the unperturbed flow. Section 
summarizes our findings. Details of linearization of the viscous 
force in the dynamical equations are given in the Appendix. 


2 DERIVATION OF THE DISPERSION EQUATION 

We start with the linear analysis of hydrodynamic equations. The 
fluid viscosity and thermal conductivity is taken into account 
through kinematic viscosity coefficient v and heat conductivity co¬ 
efficient K, respectively. 


2.1 Basic equations 

The system of hydrodynamic equations reads: 

(i) mass conservation equation 

^-hV-(p«) = 0, (1) 

Ot 

In cylindrical coordinates for axially symmetric flows: 

1 dipru,) d(pu.) 

V • ipu) = ----I- —— (2) 

r or oz 

(ii) Navier-Stokes equation including gravity force 

du 1 

— -I- (mV) •m = —V p - V<j)„ + N . (3) 

at p 

Here = —GMjr is the Newtonian gravitational potential of the 
central body with mass M, N is the viscous force. In cylindrical 
coordinates for axisymmetric flows: 

dUr dUr dUr 

- + Uf - + Ut - — — 

dt dr “ dz r 


du^ du^ du^ 


du^ du- du- 

+ Ml- 

dt dr dz 

The linearized viscous force components are specified in Appendix 
A. 

(iii) energy equation 
pHT 



(4) 

dr p dr 

- = ^d) ^ 

(5) 

r 


z p dz 

(6) 


ds 

— + (mV) • j 
dt 


= gv. 


V-F. 


(7) 


where i is the specific entropy per particle, gvisc is the viscous dis¬ 
sipation rate per unit volume, K is the universal gas constant, p is 
the molecular weight, T is the temperature, and terms on the right 
stand for the viscous energy production and the heat conductivity 
energy flux F, respectively. The energy flux due to the heat con¬ 
ductivity is 

V • F = V(-/fVr) = -K^T -Wk-WT. (8) 

Note that both electrons and photons, and at low temperatures neu¬ 
tral atoms, can contribute to the heat conductivity (see Section 
below). 

(iv) equation of state 

The equation of state for a perfect gas is convenient to write in 
the form: 

p = , (9) 

where F is a constant, cv - I /(y-1) is the specific volume heat ca¬ 
pacity and 7 = Cpjcv is the adiabatic index (5/3 for the monoatomic 
gas). 


2.2 Linearization of basic equations 

We will consider small axially symmetric perturbations in the 
WKB approximation with space-time dependence 
where r, z, 0 are cylindrical coordinates. The velocity perturbations 
are u = {Ur,u^,u^). The density, pressure, temperature and en¬ 
tropy perturbations are pi, pi, T[, and j, over the unperturbed 
values po, po. To, and ^o, respectively. As a simplification, to fil¬ 
ter out acoustic oscillations arising from the restoring pressure 
force, we will use the Boussinesq approximation, i.e. consider 
incompressible gas motion V • m = 0. In the energy equation 
we will neglect Eulerian pressure variations, pi{t,r,(j),z) = 0 
(see the justification below), but Lagrangian pressure variations 
5p{t,r(to),^{to,z{to)) are non-zero. (We remind that for infinitesi¬ 
mally small shifts a perturbed gas parcel acquires the pressu r e equa l 
to that of the ambient med i um; s ee e.e. lSniegel & VeronisI (Il960h : 
iKundu. Cohen & Dowlind (12012h for discussion of the Boussinesq 
approximation). We stress that we investigate the motion of ax¬ 
isymmetric transverse perturbations, i.e. small perturbations in the 
direction normal to the wave vector k. 


2.2.1 Dynamical equations 

In the linear approximation, the system of differential hydrody¬ 
namic equations is reduced to the following system of algebraic 
equations. 

a) . The Boussinesq approximation for gas velocity m is V • m = 

0 : 

krUr + k-U, - 0 . (10) 

b) . The radial, azimuthal and vertical components of the 
Navier-Stokes momentum equation are, respectively: 

iiOUr - 2Q.US = ikr— - - vPurlR] , (11) 

Po Po dr 

where the factor [F] takes into account the dependence of the vis¬ 
cosity coefficient on temperature q ~ F“™c = 5/2 for fully 

ionized gas and a,,,sc = 1/2 for neutral gas) in the perturbed vis¬ 
cous force component A/- (see Eq. l lA16t in Appendix A); 

ioju,!, + = -vPu^plfb ], ( 12 ) 
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where the factor [®] takes into account variations of the viscos¬ 
ity coefficient in the perturbed viscous force component (see 
Eq. llA21b in Appendix A); 




Py _ Pi dPa 
>0 pI dz 


vVuyz] 


(13) 


where the factor [Z] takes into account variations of the viscos¬ 
ity coefficient in the perturbed viscous force component N, (see 
Eq. (??) in Appendix A). Here V = kj + kj and 


= AST + r- 


dOr 

dr 


1 

dr 


(14) 


is the epicyclic frequency. For the power-law rotation ~ 
the epicyclic frequency is simply = 4 - n. In deriving these 

equations we neglected terms ~ (Kir), jkyjr) compared to terms 
~ V, see also the discussion in lAches^ ( Il978l) . 


2.2.2 Energy equation 


To specify density perturbations pi/po, the energy equation should 
be used. In the general case by varying the equation of state Eq. 
we obtain for entropy perturbations: 


^ - £l + 81 

Po Cv Po 


(15) 


On the other hand, from the equation of state for ideal gas in the 
form p = pKTIp, we find for small temperature perturbations we 
have: 


El - Pi + I 1 

Po Po To 


(16) 


Substitution of Eq. 03 and Eq. 03 into equations of mo¬ 
tion Eq. OD and Eq. 03 immediately shows that the terms with 
(pi/po)kr and (pilpo)ky in the dynamical equations are larger by 
factors rk,. and rk, than terms with pi/po arisen from the energy 
equation. This means that in the energy equation we can set Eule- 
rian pressure perturbations equal to zero, as is usually assumed in 
the Boussinesq approximation, i.e. 


-h 7 - = 0 . 

Cv Po 


(17) 


Py ^ _'Ty 
Po Tq 


(18) 


We repeat again that the Eulerian pressure variations should be re¬ 
tained in the equations of motion 03-03- Eq. 03 implies that in 
the axially symmetric waves considered here the density variations 
are in counter-phase with temperature variations. 

The viscous dissipative function 2visc [erg cm“^ s“'] can be 
written as Qvisc = pv®, where the function ® in polar coordinates 
is 


® : 


( 7 )- 


+ r 


1 dUr 






(19) 


All terms but one in this function are quadratic in small velocity 
perturbations; this term has the form: 


vp 




( 20 ) 


Writing for the azimuthal velocity u,p = u,pfl + (here for the pur¬ 
poses of this paragraph and only here we specially mark the unper¬ 
turbed velocity with index 0, not to be confused with our notations 
for perturbed velocity in Eq. <1 Ib -Eg. l ll2b above and below), we 
obtain for the viscous dissipation function 


dO. 

Qvhc = vpr-r- 
dr 


r— - zikyU^i 

dr 



-I-quadratic terms. (21) 


Here fl = u^^/r is the angular (Keplerian) velocity of the unper¬ 
turbed flow. The first term in parentheses describes the viscous en¬ 
ergy release in the unperturbed Keplerian flow. For this unperturbed 
flow we have 




[ridn/dr)]^ 

Wo 


9 

” A^'^W 


( 22 ) 


Thus, the entropy of the unperturbed flow changes along the radius. 
However, on the scale of the order of or smaller than the disc thick¬ 
ness zo. the entropy gradient can be neglected. The second term in 
Eq. mj vanishes if = 0 (and then there are no viscous dissipa¬ 
tion effects to linear order), therefore we will consider only two- 
dimensional transverse perturbations with k, t 0,ky t 0. 

We emphasize that in our analysis we neglect the background 
entropy gradients, which can be present in real flows, i.e. we will 
consider the flow in the local neutral equilibrium. This is done to 
exclude the effects of these gradients on the evolution of small per¬ 
turbations. As is well known, with inclusion of the background 
entropy gradients, the Brunt-Vaisala frequencies arise. If their 
squares are positive, they stabilize perturbations. If their squares 
are negative, they signal emergence o f convection (see, for exam- 
ple. lKato. Fukue & Mineshiga l ll99^ . for more detail). 


The right side of the heat conductivity equation Eq. l[8j for 
small temperature perturbations, with account for the dependence 
of the thermal conductivity coefficient on temperature and density 
K ~ T‘^p‘‘, in the linear order can be recast to the form 


kAT + Wk-WT = -kE^To^ - 

(c - d)KTo ^ iky ^ ^ - (c - d)KTo ^ iky ^ fi, (23) 

(here we have used the relation l llSb ). It is useful to rewrite the 
right side of this equation in the form 


To 


kr 1 ^T^ 


ky 1 dTo 


- kVTo—\ \ + ic - dli^T -^ -H(c-d)i-f-^ (24) 


V To dr 


V To dz 


to see that in so far as ~ - I/r, - l/zo and that terms 

Tq dr ' ’ Tq dz ' ^ 

~ k,./r should be neglected compared to terms ~ V, only the first 
term ~ k^ and third term k,/zo should be retained in this equation. 
Therefore, the energy equation Eq. 0 turns into 

. PoKTo dO. ,2rrT\yy,y 

icj - =—likyvpor - Us — Kk lo — [cj, (2h) 

p dr To 

where we have introduced the correction factor 


[£] = 1 -H ;(c - d) 



To dz ■ 


Like in the linearized continuity equation V • « = 0, here we 
have neglected the term u,plr. The first term in the right side of 
Eq. l l25b corresponds to the energy generation in axially symmet¬ 
ric sheared flow due to viscosity, and the second term means the 
entropy perturbation smoothing due to heat conductivity. The first 
term (ky/rliu^oluy^q, and the second term Vq/Pr, where the 
Prandtl number is the ratio of the heat conductivity to dynamical 
viscosity coefficient. While the k^/r is small relative to V, the coef¬ 
ficient (u^fllu^y is very large for thin discs, and therefore the heat 
generation term should be retained in the energy equation. 
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2.3 Dispersion equation 

By substituting Eq. il7t and Eq. l ll 8 t into Eq. i25b . we find the 
relation between the density variations and in the Boussinesq 
limit with zero background entropy gradients: 


Po 


kP[E] \ 2ikrVr(dQ.Idr) 
CpPoKlij) CpKTolfi 


(26) 


Here Cp = ycy = 7 /( 7 -1) is the specific heat capacity (per particle) 
at constant pressure. 

It is convenient to introduce the dimensionless Prandtl num¬ 
ber: 


Pr = 


vpaCp vpo(1ilp)Cp vpoCRIp) 7 


7 - 1 


(27) 


The Prandtl number defined by Eq. J27b for fully ionized hydrogen 
gas (7 = 5/3), where the heat conductio n is determined by light 
electrons, is quite low fsee lSDit^ ( ll962l) ): 


Pte 


0.406 


20 • 0.4 • 0.225 • (2/jr)3/2 


1/2 


: 0.052. 


(28) 


Note also that in this case <7 = 0 and c = 5/2 in the heat conductivity 
coefficient. 

By expressing the heat conductivity coefficient k through kine¬ 
matic viscosity coefficient v using Eq. ( 127b and after substituting 
Eq. l l26b into Eq. dl lb . we arrive at: 


P\ 2 ^ 

ikr— = (tOJ + vk [R\)Ur -I- --Ur 

Po {i(jj + vk^[<P\) 

1 1 (9po x^ivkr(d\nQ.ld\Tir) 

Cp Po dr (io) -I- vF[il)])(iw -1- vk?-[E\l'Pr) 


(29) 


Now by substituting Eq. ( l26t into Eq. with account for 
Eq. dlOb . we arrive at 

ik,— = — -{iu) + vfc^[Z])Ur 
' po k. 


1 1 dpo K^ivky{d\nQ.ld\nr) 

Cp Po dz (iai -(- vF[0])(iw - 1 - v/:2[£]/Pr)^ 


(30) 


Finally, by subtracting Eq. d30b multiplied by k,. from Eq. d29b 
multiplied by k- we arrive at the dispersion equation: 


(ioj + 


(idj + vP[R])-j^ + (iw + 



2 

llj 



7 - 1 


ikr 


where 


7 (io) + vkHE]IPi) 
c/lnfl 1 dpo 


A-^B 
k. 


A = V- 


B = v 


dlnr Po dr 
<71nfi 1 dpo 


0, (31) 


(32) 


(33) 


dlnr Po dz 

The expression in the square brackets in Eq. (EB above can be 
rewritten in the equivalent form: 


1 + 


7 - 1 


dXnElldlnr 


, _Pr 

^rSoeff ^ Sz 


7 (iw + yk2[£]/Pr) KTo/p 

(34) 

where gr.^// = -llPoidpoldr) and g, = -llpo(dpoldz) are the ef¬ 
fective radial and vertical gravity accelerations in the unperturbed 
flow, respectively. Clearly, the term Kg^^ff ~ k^/r is much smaller 
than {kf/k^)g^ ~ {Pjk^XIzo and will be neglected in the further 


analysis. Note that if the dynamic viscosity coefficient is indepen¬ 
dent of temperature (i.e. = 0 ), correction factors [7?] = [O] = 

[Z] = 1 and the first line in Eq. d31b is simplified to {icj + vP)^. We 
will see below that deviations of these correction factors from unity 
insignificantly affect the result of our analysis. 


3 ANALYSIS OF THE DISPERSION EQUATION 

First consider the limiting case where A = B = 0 and [7?] = [O] = 
[Z] = 1, i.e. the case where the viscous energy generation in the 
energy equation is ignored, but the viscosity is retained in the equa¬ 
tions of motion. Then Eq. d31b represents the well-known Rayleigh 
dispersion equation for viscous fluid: 

(itu + vk^)^ = - (35) 

which describes two Rayleigh modes. Depending on the sign of 
these modes are oscillating (if the epicyclic frequency > 0 , the 
angular momentum increases with radius), or exponentially grow¬ 
ing (the unstable Rayleigh mode) and exponentially decaying (if 
the epicyclic x^ < 0 , the angular momentum decreases with ra¬ 
dius). The arising of the unstable Rayleigh mode corresponds to 
the classical Rayleigh criterion of instability of a shearing flow. As 
can be easily seen from Eq. l l35b . the viscosity stabilizes the unsta¬ 
ble mode at short wavelengths (large k). Two solutions of Eq. ([35} 
for typical viscosity parameters discussed below and (the 

Keplerian motion) are shown in the fourth column of Fig.[T] 

However, in the general viscous case where A,B t 0, the dis¬ 
persion equation Eq. I l31b turns into a cubic equation, that is, the 
third mode arises (the X-mode) due to the viscous heating of the 
fluid. At non-zero A and B, one of the Rayleigh modes (which is 
stable in the dissipationless case) becomes exponentially unstable 
in a wide range of wavenumbers. We stress that these modes remain 
stable for either kr = 0 or k- = 0. Indeed, the dispersion equation for 
perturbations with = 0 is reduced to Eq. ( 1351 above. Eor pertur¬ 
bations with k. = 0 , the dispersion equation turns into ioj + vkj = 0 , 
i.e. is reduced to an exponentially decaying standing wave. We 
stress that in our case the larger viscosity, the higher instability in¬ 
crement. This is opposite to the situation where a poloidal magnetic 
field is present, when the magneto-rotational instability is devel¬ 
oped: increasing viscosity decreases the MRI increment. Clearly, 
there is no viscous instability of the Rayleigh modes in the invis- 
cid case (v = 0 ) or in the shearless case (solid-body rotation with 
n = 0 ). 


3.1 Dimensionless dispersion equation 


For numerical analysis, the cubic dispersion equation Eq. 1(33 can 
be conveniently rewritten in the dimensionless form. To do this, 
we multiply Eq. l l31l l through the factor (iu) + vk^[7i]/Pr), divide 
the obtained equation through Q? and introduce new dimensionless 
variables: 


a> 



kr. 


,2 _ X 

' " f?- 


(36) 


The kinematic viscosity is v = lu„ where I is the effective mean 
free path of ions, m, = ojyRTolp is the characteristic velocity in 
the unperturbed flow which is about thermal velocity of ions, so 
the dimensionless combination vPiQ. becomes: 


vP 

IT 


= a(kr)^ . 


(37) 
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n=3, Pr =0.052, u /u =0.01, l/r=0.01, k /k=0.999 
e s 0 r 



ii=3,u^/y^=O.O1,l/r=0.OU/k=O.999,UO 



n=3,y/y40Ur40UM9 




R+ 



R- 




0 5 10 15 20 25 30 

kr 





Figure 1. Three modes of dispersion equation for electron thermal conductivity in fully ionized plasma (the Prandtl number Pr£.=0.052). R+ and R- mark 
two Rayleigh modes, one of which (R+ if > 0 and R- if < 0) becomes viscously unstable. Upper row: Re^j, bottom row: Im^. The first two column 
show the case krik = 0.999,k- > 0. The left column zooms the region of small kr to see the behaviour of three modes as kr —> 0: Re and Im parts of the 
X-mode always starts from zero, while at = 0 Re and Im parts of the Rayleigh modes is non-zero and zero, respectively. The third column shows the case 
krIk = 0.999, k- < 0: the imaginaiy part of the modes is unchanged and the real part changes sign. For comparison, the fourth column shows two (stable) 
Rayleigh modes as the solution of Eq. (35). 


Here we have introduced the dimensionless coefficient 


Here the dimensionless factors [iJ], [®], [Z] and [£] have the form 


a = 



(38) 


Formally, 1/a is the Reynolds number defined as Re= (u^r)lv, but 
as we will see below, for a specified Reynolds number, different 
solutions are realized. 

The vertical pressure gradient in coefficient B in Eq. OH turns 

into 


1 dpo _1_ 

Po dz Zo ’ 


(39) 


where zo is the characteristic disc height. Using the relation for thin 
accretion discs 


- = VifTr 

r 



(40) 


(where the dimensionless coe fficient Hi takes into acco unt the 
model vertical disc structure, see lKetsaris & Shakural il998l) : in nu¬ 
merical calculation below we shall assume -^Hi/y = 2), we obtain 
the dispersion equation in the dimensionless form: 


(iS) + a{krf-[^Y) 


1^2 

(iai + a(fcr)^[R])-| -l- (io) + a(kr)~[Z]) 




'k^ 


T U: 


l\ k, 


1 -'t( 7-1) Vni/r ,, x2rcn/T3 
2 10) + a{kry[E\IV\: 


= 0. (41) 


[Z] = [l-i 2 a„..(|)^(^)], 

[£] = [l-i(c-rf)(|)^(^)]. (42) 

The inspection of Eq. j41b reveals the following properties of 
the solution: 

• the solution should be independent on the radial direction of 
the perturbation wave since the radial component of the wave vec¬ 
tor appears as Change of the sign of k. reverses the sign of the 
real part of the solutions (see the second and third column in Fig. 

0 ; 

• in the limit of small viscosity, the second bracket in Eq. (EB 
becomes real in the first order, suggesting the stability. The de¬ 
crease in the viscous instability increment with decreasing Ijr is 
clearly seen in Fig.[3j 

• for uju^ ~ 0.01 (thin discs) and small k-/k « 1, k,- ~ 1 
and kr ~ 10, where the viscous instability appears (see below), the 
most appreciable correction is for the [R]-factor. However, in the 
dispersion equation ED the term (iw -l- a{kr)^[R]) is multiplied by 
the small value (k-lk)^, and therefore the effects from the correction 
factors [R] - [£] on the solution of the dispersion equation should 
be not significant, as indeed we found to be the case. 
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This dimensionless dispersion equation for w as a function 
of the dimensionless wavenumber {kr) is to be solved for differ¬ 
ent values of the dimensionless parameters: the Prandtl number 
Pr, which characterizes the effect of thermal conductivity, Ijr and 
u ju^, which describe the viscosity, and k^lk, which determines the 
direction of the wave (evidently, (k^lk)^ = 1 - (kr/kf). 


3.2 Numerical solution of the dispersion equation 


In principle, it is possible to carry out analytical investigation of 
the properties of the solutions of the cubic equation Eg. 
in a way similar to study of MRI modes bv lPessah & Chan 
However, the main aim of the present paper is to show the ex¬ 
istence of the viscous instability in shearing flows, Iherefore we 
will numerically solve Eq. (ED for different representative param¬ 
eters. Everywhere below in this Section we shall consider the phe¬ 
nomenologically important Keplerian case with n = 3 and % = 1. 
This does not restrict our analysis, since the instability persists at 
any n but n = 0 (see the next Section). 


3.2.1 Case of electron heat conductivity 

We start with the electron heat conductivity in a fully ion¬ 
ized plasma. We remind that in this case the Prandtl number is 
PiV=0.052, the dynamical viscosity coefficient is ~ the heat 
conductivity coefficient isK ~ so that aviic = 5/2, c = 5/2 and 
d = 0 in Eq. ED. Fig. [D shows the real (upper panels) and imagi¬ 
nary (bottom panels) parts of three solutions of the cubic dispersion 
equation Eq. ED w as a function of the dimensionless wavenum¬ 
ber kr. All three solutions of this equation are complex since the 
dispersion equation has complex coefficients. The dimensionless 
parameters are uju^ = 0.01 (thin discs), l/r = 0.01 (the maximum 
possible free-path length of ions, not to exceed the disc thickness), 
krik = 0.999 (the direction of perturbations with increment close 
to maximal one for these parameters, see Fig. |D- Two Rayleigh 
modes modified by viscosity are marked as R+ and R-, according 
to the sign of their real parts at kr —> 0. The first two columns show 
the solutions for kr = 0.999 and positive k- > 0 and negative k, < 0, 
respectively. It is seen that the sign of k- determines which of the 
Rayleigh modes, R-l- or R-, becomes unstable. It is also seen the un¬ 
stable mode is that which has the real part intersecting with the new 
X-mode (the latter is always stable, i.e. has a non-negative imagi¬ 
nary part, representing an oscillating wave). The unstable Rayleigh 
mode has a non-zero increment already for long perturbations with 
kr —> 0. It has a maximum increment of ~ 0.1 at fcr sj 5 and is 
stabilized by viscosity for kr >25. 

Fig.[2]illustrates the effect of changing the perturbation prop¬ 
agation wavevector value kr/k in the range from 0.9 to 0.9999. 
It is seen that at kr/k = 0.999 the instability increment is about 
maximum (we did not investigate the exact value of krjk for maxi¬ 
mum increment, which, if necessary, can be straightforwardly done 
by differentiating the dispersion equation with respect to kr/k and 
equating the result to zero). 

Fig. [D shows the unstable Rayleigh R+ mode behaviour with 
changing the viscosity parameter Ijr and other parameters fixed as 
in Fig. ID It is seen that diminishing the particle free-path length 
from the maximum possible value (Ijr = 0.01 in this case) by an 
order of magnitude decreases the Rayleigh R+ mode instability in¬ 
crement by about two times, but increases the instability interval 
from kr 25 to kr == 70. 

Fig.lDalso illustrates the effect of increasing or decreasing the 


n=3,u /u =0.01, l/r=0.01 

S (j) 



Figure 2. Imaginary pai't of the viscously unstable mode R+ in fully ionized 
gas with electron heat conductivity (Pre=0.052) and viscosity parameters 
= 0.01,//r = 0.01 for four values of the wave vector k,-/k = 0.9, 
0.99, 0.999, and 0.9999. 


disc thickness by three times {uju^ - 0.03 and uju^ - 0.003, re¬ 
spectively). The imaginary part of the unstable R-l- mode is shown 
for three values of the viscosity parameter Ijr = 0.03, 0.01 and 
0.003. Note the close similarity (almost identity) of the curves 
to those in Fig. [D but the stretching of the kr variable by three 
times. This reflects an almost self-similarity of the dispersion equa¬ 
tion Eq. ED with respect to the dimensionless viscosity a{krf = 
(uju^){llr){krf. 


3.2.2 Case of radiative heat conductivity 

For a mixture of electrons and photons, the heat conductivity can 
be characterized of an effective Prandtl number defined as 


Pr Pr,, Pr Pr 


Pr 


- = — + - = -! + — = —l + - 


Pr. 


Pr. 


qy 


where the heat flux due fo elecfrons is 
1 

and fhe heat flux due to photons is 

qy = -^clyW{aJ‘*) 

where a^ is the radiation constant. Therefore, 

‘ly „ c Zo 
— -(3 - 

qe Us T 


(43) 


(44) 


(45) 


( 46 ) 
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n=3, u^/u^=0.03, kyk=0.999 n=3, u^/u^=0.01, kyk=0.999 n=3,u^/u^=9.003, k/k=0.999 



Figure 3. Almost precise self-similarity of the solution for different disc thickness parameters Shown are three cases for Pr5=0.052, kyjk = 0.999 and 

(from left to right) u juphi = 0.03,0.01,0.003 with different values of l/r. It is seen that the combination a{kr)^ is almost exactly the same for all three figures. 


where ft = Pyjpg^ is the radiation to gas pressure ratio, t is the 
effective optical thickness of the disc and cTj, is the electron-ion 
interaction cross-section. Noticing that n^cTgjU,, = v^,- = is 

the electron-ion collisional frequency, Eq. gill can he recast into 
the form 

‘ly ^ fl_ c 

2jt T Us 

Fig.glshows the effect of decreasing the effective Prandtl num¬ 
ber due to increase of the radiation heat conductivity. Two cases 
with Pr=Pr5/2, Pr^/ll are shown in comparison with the case of 
electron heat conductivity only. It is seen that the radiation con¬ 
ductivity in fully ionized plasma strongly decreases (hut does not 
vanish) the instability increment. 

3.2.3 Case of cold neutral gas 

Let us discuss the case of cold neutral gas. In this case the 
Prandtl number Pr„=2/3 according to simplified kinetic theory 


jHirschfelder. Curtiss & Birdll954 

) and the heat conductivity 

efficient k ~ 7^2 (c = 1/2, rf = 0) ( 

Spitzerll 19621). 


Fig.m shows the imaginary part of the unstable mode R-l- for 
the standard parameters = 0.01, krjk = 0.999 used above 
in the case of ideal neutral hydrogen gas with Pr„ = 2/3 and the 
viscosity parameters l/r = 0.01 and l/r = 0.03. As above, the de¬ 
crease in the particle free-path length widens the instability wave¬ 
length interval and decreases the instability increment. In this case, 
the instability increment is maximum at kr == 17 and is about 0.18, 
almost two times as large as in the case of the purely electron heat 
conductivity in fully ionized gas discussed above. Therefore, the 
viscous instability turns out to be the most strong in the case of 
cold neutral gases. 

4 SHEARED FLOWS WITH NON-KEPLERIAN 
ROTATION 

Here we discuss the behaviour of the viscously unstable Rayleigh 
mode in flows with possible non-Keplerian rotation (i.e where oc 


r^" and n 3). The solid-body rotation case with n = 0 was already 
discussed above. In that case there is no shear and the coefficients 
A = S = 0 in Eq. Oil but the viscosity remains in equations of 
motion (see the discussion at the beginning of Section^. 

The case of a Rayleigh-unstable flow with n > 4 (i.e. with 
specific angular momentum decreasing outward) is shown in Fig. 
1^ Here the R- mode is unstable (unlike in the Keplerian case with 
k- > 0) in a wide range of kr with non-zero negative imaginary part 
at fcr —» 0 and the maximum instability increment ~ 0.2. 

Now consider a flow with increasing angular velocity with ra- 
dius, i.e. with n < 0. As is well know n, such flows are MRI-stable 
dVelikhovH 19591 : Ichandrasekhailll96(il) . However, the viscous insta¬ 
bility discussed in this paper persists in this case (see Fig.|7). Like 
in the case with n = 6, the R- mode is unstable in a wide range of 
kr. 

Finally, the flow with constant angular momentum {n = 4) 
corresponds to x = 0 and deserves special consideration. Such 
flows can be realized in various astrophysical situations, e.g. in 
quasi -spherical accretion with angular momentum onto compact 
stars dShakura et alj|20f^ . If the correction factors [R] - [£] I I42I 
were ignored, pure decay of perturbations due to viscosity would 
take place, w = ivV. However, if they are taken into account, the 
solution of Eq. | |4H is 

a> = ivV + OsiscvV Y 77^ — > (48) 

k {krf Us 

representing decaying oscillations with frequency ~ asiscCliflkf), 
which can be much smaller than the Keplerian one. 


5 DISCUSSION 

5.1 Justification of the approximation of incompressibility 

As is well known (see lLandau & Lifshitzl ( IT95^ ). the approxima¬ 
tion of incompressibility requires the characteristic time of the den¬ 
sity change in a fluid to satisfy the relation t » L/Cj, where L is 
the characteristic scale of the problem. For perturbations with the 
characteristic frequency o) = Inlr and wavenumber k = 2jr/L this 
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n=3,u /u =0.01, k/k=0.999 



Figure 4. Imaginary part of the unstable mode R+ for = 0.01, //r = 
0.01, krik = 0.999 and dilferent values of the effective Prandtl number ED. 
illustrating the effect of radiative conductivity growth. 


n=3,Pr =2/3, u /u =0.01, k/k=0.999 



Figure 5. Imaginary part of the unstable mode R+ in the case of ideal neu¬ 
tral hydrogen gas with Pr„ = 2/3 for krIk = 0.999, Usju^ = 0.01, and 
viscosity parameter Ijr = 0.01. 


general relation yields o) and in the thin discs with ~ Hzo 

we obtain the condition of the incompressibility in the form 

feo » ^ • (49) 

Writing feo = (kr)(zolr) = (kr)(uju^), this condition becomes 


(kr) » 


(cdIQ) 

(uju^) 


(50) 


For thin discs with uju^ ~ 0.01 - 0.03 and for the found mode 
frequencies w < O.lfl we see that the assumption of the incom¬ 
pressibility is valid for modes with {kr) » 3 - 10. This implies that 
in the range (kr) ~ 10-50 where the viscous instability considered 
here reaches maximum increments (especially in the case of cold 
neutral gases) the assumption of incompressibility is justified and 
sound wave modes can be ignored. 


5.2 Damping by entropy gradients 

So far we have ignored the possible radial and vertical entropy gra¬ 
dients, i.e . have dealt with locally ad i abatic flow. As is well known 
(see, e.g.. lKato. Fukue & Mineshigd l ll998[) ). the presence of non¬ 
zero entropy gradients s, = dsjdr and = ds/dz can stabilize 
instabilities. For example, if the vertical temperature gradient in 
a flow is non-adiabatic, dT/dz < dTjdzad = 8zlCp, the restor¬ 
ing gravity force would suppress the development of convection, 
leading to an oscillatory vertical motion of a gas parcel with the 


Brunt-Vaisala frequency N^. Qualitatively, it is expected that if this 
frequency is larger than the instability increment, the perturbation 
amplitude will not increase. To quantify this, we introduce the en¬ 
tropy gradients into the right-hand side of energy equation Q, and 
arrive at the modified dispersion equation: 


(ioj + (ioj + vlc^[R])i^ + {ioj -t- y/:^[Z])^ 


k. 


1 -I- 


(iw -k 




(icj + yk2[£']/Pr) 



ikr 

(iw + yk2[£]/Pr) 



: 0 . 


(51) 


Here = -Srgr and = -S .g, are the Brunt-Vaisala frequen¬ 
cies. It is seen that it is the vertical Brunt-Vaisala frequency that 
mostly affects the results, the radial oscillations being suppressed 
by small factor k^/k^. We find that in the case of Keplerian rotation 
of ionized ideal gas with Pi>=0.052 and k^/k = 0.999 the viscous 
instability discussed above disappears for N^jQ. > 0.3. Neutral gas 
with Pr„=2/3 is stabilized if > 0.35. For example, for a poly¬ 
tropic thin accretion discs with vertical struct ure described by the 
polytr opic index n' , P = , discussed in iKetsaris & Shakural 

(Il998l) . NI = 2z^{n' — 3/2)n|/(l - z^), where is the Keplerian 
rotation frequency. The Brunt-Vaisala frequency N, averaged over 
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n=6, u /u =0.01, l/r=0.01, k/k=0.999 



n=-2, u /u =0.01, l/r=0.01, k/k=0.999 

s $ r 





Figure 6. Three modes of dispersion equation Eq. <4U in the case of a 
Rayleigh-unstable flow with n = 6. 


Figure 7. Three modes of dispersion equation Eq. | |4U in the case of a flow 
with angulai' velocity linearly increasing with radius 0. ~ r {n = -2), which 
is stable against MRI. 


the disc height zo is (NJCIk) = V(”' “ 3/2)/2. Therefore, the value 
N, = 0.3 corresponds to a polytropic index n' == 1.7. Of course, 
realistic flows can he not polytropic, and therefore effects of the 
entropy gradients on the viscous instability should be investigated 
separately in each particular case. 


6 SUMMARY AND CONCLUSION 

In the present paper we have performed a linear local WKB analy¬ 
sis of time evolution of small axisymmetric perturbations in sheared 
hydrodynamic laminar flows. As a simplification, we have used the 
Boussinesq approximation for the description of the perturbations, 
but included the viscous dissipation and heat conductivity terms in 
the energy equation. This procedure led us to a third-order algebraic 
dispersion equation with complex coefficients (see Eq. I I41U . The 
inclusion of these terms makes one of the Rayleigh modes (with 
positive or negative real part depending on the sign of the wavevec- 
tor component U) unstable for long-wave perturbations for locally 
adiabatic case (i.e. ignoring local entropy gradients). The new X- 
mode of this cubic equation is found to be always stable (i.e. has a 
positive imaginary part). 

We have studied numerically the behaviour of the unsta¬ 
ble Rayleigh mode in the most interesting case of thin Keple- 
rian accretion discs for different values of the viscosity (which is 
parametrized by the mean free-path length of ions), disc thickness 
(which is described by the ratio of the sound velocity to the unper¬ 


turbed tangential velocity in the flow), the directions of the pertur¬ 
bation propagation (which is described by the ratio of wave vector 
components kr/k-), and the Prandtl numbers (which describe the 
heat conductivity effects). We have found that the value of heat 
conductivity mostly affect the instability increment, which is found 
to be maximum 0.2 of the local Keplerian frequency in the case 
of cold neutral gas with the highest value of the Prandtl number 
Pr„ = 2/3 (see Fig.j^. In the fully ionized gas characterized by the 
Prandtl number Pr,, = 0.052 for purely electron heat conductivity, 
the instability increment is about 0.1 and decreases with increasing 
the role of the radiation heat conductivity (Fig.|4ll. The instability 
increment does not sensitive to the direction of propagation of per¬ 
turbations (the sign of wavenumbers kr and kj (Fig.[T) and persists 
as long as shear and viscosity are present in the flow and the flow is 
not iso-momentum when the epicyclic frequency vanishes, i.e. for 
any law of the angular momentum iV ~ r^" (see Fig.|^and Fig.|7]l. 

In the presence of viscous dissipation, the instability arises 
when the pressure gradients along radial or vertical coordinates are 
non-zero, suggesting its convective nature: the heat generation in a 
sheared viscous flow in the gravity field of the central star makes the 
flow convectively unstable. Different aspects of convection in cold 
accretion discs, especially suitable for the physics of protoplanetary 
discs, has be en addressed in many pa pers, start ing from the pio¬ 
neer p a per bv | Lin & Pap aloizoul l ll980ti (see also iRvu & Goodman! 
l ll99i) : lLesur & Qgilviel ( l2oiT and references therein). 

We show that the incompressibility approximation is applica- 
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ble to describe small perturbations in thin accretion discs with not 
very long wavelength (kr » 3 - 10). At longer wavelengths, acous¬ 
tic perturbations should be taken into account. On the other hand, 
the local WKB analysis is applicable only for kr » 1. Thus, the 
found instability with maximum increment at kr ~ 10-50 seems 
to be robust under our assumptions. 

Thus we conclude that the viscous instability of one of the 
classical Rayleigh mode discovered in the present paper may he a 
seed for the development of turbulence in sheared flows which are 
hydrodynamically stable according to the classical Rayleigh crite¬ 
rion (i.e. in which the angular momentum increases with radius), or 
stable against MRI (e.g. flows with angular velocity increasing with 
radius). This instability is certainly worth investigating further. 
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APPENDIX A: LINEARIZATION OF VISCOUS FORCE IN DYNAMICAL EQUATIONS 

In cylindrical coordinates for axisymmetric flows the viscous force components read (see e.g. lKato. Fukue & Mineshi^ ( Il998ll l 

Kr _ ^ 

p\r dr r dz 


(Al) 


^ / _l_ djr^tr^) ^ dCjp 


p \r 


dr 


dz 


(A2) 


The viscous stress tensor components are: 


At = 1 / 1 djrtr^) 

^ p\r dr dz 


duy 2 


tr<p = pr 


d(u^lr) 

dr 


du, dUr 


(A3) 

(A4) 

(A5) 

(A6) 


Ur 2 

= 27/— + (f - -77)V • u , 
r 3 


= '/■ 


dUrf, 

dz 


du, 2 

t^z = 2;/-^ +{(- -77)V • « . 

(Here ij = pvis the dynamical viscosity, ^ is the second viscosity.). 

Below unperturbed and perturbed components will be marked with indexes 0 and 1, respectively, and therefore 


hi'f Uf 2 ^ 


P-Po+Pli 7/— 7/0 +7/1 > T—Tq + Ti. 


For perturbed variables pi, uir, 4 ,,z),i’ Ty taken in the form of plane waves ~ exp(im? - k,./- - Cz) the partial derivatives simply 
djdr = —ikr, didz = —ik.. The dynamic viscosity of interest here is a function of temperature only, rj ~ therefore 

dp Ti 770 dTo 

o “ ^viscPO rj, ^kr + (Trisc „ , 

dr To To dr 

dp _ 2’i po dTo 

r, ~ ^ViScPO rj, + O^visc ™ 5 

dz To To dz 

(By varying 77, we neglected logarithmic dependence on temperature and density in the Coulomb logarithm). 


(A7) 

(A8) 

(A9) 

(A 10) 
becomes 

(All) 

(A12) 


Al Radial component 

Inserting the stress tensor components into Eq. JaB yields 


1 

Nr=- 

P 


cfUr d^Ur 77 dUr 


dp dur dp du- dp duy 


^ dr- dz^ r dr ^ r^ dr dr dz dr dz dz 


After linearizing we find: 


K.i = - 
Po 


,2 T, PodTo PodTo PodTo 

io or To dz To dz 


(A13) 


(A14) 


The second term is ~ and is small compared to the first term ~ V and last two terms ~ k^/z (we remind that in thin discs considered here 
Zo/r ~ uju^pfi « 1). Noticing that + Ur^ikr = 0 from the continuity equation, we obtain: 

Nr,i=-vVur,i[R\ (A15) 

where the factor that takes into account the dependence of viscosity on temperature is defined as 


V -V 

[P] = 1 + i^-r-rfl-Ur, 


_L^ 

'To dz 


(A16) 
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A2 Tangential component 


Substituting stress tensor components into Eq. IA2t yields: 



rj du^ drj du^ 

^ dn ^ dz^ r dr ^ r^ dr dr 


drj drj du^ 
dr r dz dz 


The linearizing leads to: 


A/’ai = - 
Po 


Ti rjoTi\du^,o 

^visc^O ^ ^visc ^ 

To r To 


dr 


rjoTi 2 ■, 1 ^To 

- O^i'isc^ T^r^ip.O — t7o^“(6,l “ “-770—i 

r^ To To dz 


-J_£l 

Po Po 


Po dTp ijo \ du^,o l oyisc PodTo po\ 

T„ ar r ] ar ^ \ r T„ ar r2 


(A17) 


(A18) 


All terms in the second square brackets are ~ 1/r^ and can be neglected compared to terms in the first square brackets. The latter can he 
rewritten as the sum of two terms: 


Kf - 1 

- —zf 
Po To 


I Po \ du^,o Po 

y ^visc^O^^r "I" ^visc j ^ ^visc 2 


, 2 “i^,i /, 1 dTo\ 

- vk^ - 1 + u^,0 ■ 

U^.o \ To dz ' 


(A19) 


Here three terms in the first brackets are ~ fcr/r, ~ 1/r^, ~ 1/r^, respectively, compared to terms ^ P in the second brackets, and hence can 
he neglected. Therefore, we are left with 


N>l>.i = [O] 


where 


[®] = 1 + 


LEl 

To dz ■ 


(A20) 


(A21) 


A3 Vertical component 


Substituting the stress tensor components into Eq. tA3l l yields: 





Ip drj\ldu. 3m,. \ 
\ r ^ 3r / \ 3r dz j 


drj du, 

+ 2 _ - _^ 

dz dz 


(A22) 


The linearization of terms in the middle brackets yields terms ^ k^jr, which are small compared to the term ~ P arisen from the second 
derivatives. Terms ~ k^jz arisen from the linearization of the last term, however, should be retained. Therefore, we finally find: 


iV,,i = -vPu,^i[Z] 


(A23) 


where 












